GAM to model the number of clicks per second as a function of depth and whale id
Carolina S. Marques
1 Reading data
#read data obtained in "01_data_processing.rmd":
load("processed_data.rda")
load("cue_rates_obtained.rda")
sampledata<-sampledata[sampledata$rsecs>=0&!is.na(sampledata$clicks),]
sampledata$time11<-as.numeric(sampledata$rsecs)
originaldata$depth<-originaldata$Depth
originaldata$time11<-as.numeric(originaldata$sst)
onedata$depth<-onedata$Depth
onedata$time11<-as.numeric(onedata$sst)2 GAM Models
2.1 GAM number of clicks per second as a function of depth by whale and as a funtion of time by whale
gamwhales1<-gam(clicks~s(depth,k=7,by=as.factor(whale))+s(time11,k=1,by=as.factor(whale)),data=sampledata,family=poisson(link="log"))## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible2.1.1 Summary of the model
summary(gamwhales1)##
## Family: poisson
## Link function: log
##
## Formula:
## clicks ~ s(depth, k = 7, by = as.factor(whale)) + s(time11, k = 1,
## by = as.factor(whale))
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.39154 0.02194 63.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(depth):as.factor(whale)Freya 4.325 4.940 17.035 0.003860 **
## s(depth):as.factor(whale)Thora 5.644 5.940 21.738 0.000932 ***
## s(depth):as.factor(whale)Mára 2.210 2.671 1.631 0.598642
## s(depth):as.factor(whale)Frida 1.000 1.000 19.735 9.49e-06 ***
## s(depth):as.factor(whale)Eistla 2.349 2.919 16.056 0.000969 ***
## s(depth):as.factor(whale)Balder 4.242 4.999 33.387 3.60e-06 ***
## s(depth):as.factor(whale)Jonas 4.750 5.447 29.933 4.59e-05 ***
## s(depth):as.factor(whale)Mutti 1.000 1.000 21.955 3.52e-06 ***
## s(time11):as.factor(whale)Freya 1.729 1.925 13.614 0.005081 **
## s(time11):as.factor(whale)Thora 1.730 1.926 4.325 0.075160 .
## s(time11):as.factor(whale)Mára 1.000 1.000 0.399 0.527811
## s(time11):as.factor(whale)Frida 1.396 1.634 24.510 9.76e-05 ***
## s(time11):as.factor(whale)Eistla 1.000 1.000 2.812 0.093580 .
## s(time11):as.factor(whale)Balder 1.000 1.000 17.759 2.56e-05 ***
## s(time11):as.factor(whale)Jonas 1.001 1.002 0.159 0.690277
## s(time11):as.factor(whale)Mutti 1.000 1.000 0.369 0.543610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0464 Deviance explained = 6.19%
## UBRE = 1.8611 Scale est. = 1 n = 12032.1.2 AIC
AIC(gamwhales1)## [1] 6732.8582.1.3 Predicting influence of the different expl. variables
dd1<-sampledata%>%
group_by(whale)%>%
summarise(depth=mean(depth,na.rm=T))# creating a data set to predict the influence of time
dd2<-sampledata%>%
group_by(whale)%>%
summarise(time11=mean(time11,na.rm=T))# creating a data set to predict the influence of depth
originaldata11<-data.frame(depth=sampledata$depth,whale=sampledata$whale)# creating a data set to predict the influence of depth
originaldata11<-full_join(originaldata11,dd2)# joining both infos to create a data set to predict the influence of depth## Joining, by = "whale"originaldata111<-data.frame(time11=sampledata$time11,whale=sampledata$whale)# creating a data set to predict the influence of time
originaldata111<-full_join(originaldata111,dd1)# joining both infos to create a data set to predict the influence of time## Joining, by = "whale"datat1<-data.frame(pclick=predict(gamwhales1,newdata = originaldata11,type="response"),depth=originaldata11$depth,whale=originaldata11$whale)# predicting the influence of depth in the number of clicks produced
confintervals<-data.frame(confintervals=predict(gamwhales1,newdata=(originaldata11),type="response",se.fit=T))# get the corresponding CI
datat1$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat1$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit
datat1$variable<-"Depth (m)"# which variable corresponds
datat11<-data.frame(pclick=predict(gamwhales1,newdata = originaldata111,type="response"),time11=originaldata111$time11,whale=originaldata111$whale)# predicting the influence of time in the number of clicks produced
confintervals<-data.frame(confintervals=predict(gamwhales1,newdata=(originaldata111),type="response",se.fit=T))# get the corresponding CI
datat11$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat11$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit
datat11$variable<-"Time"# which variable corresponds
datat<-full_join(datat1,datat11)#Joinning both infos about the variables ## Joining, by = c("pclick", "whale", "bigfit", "smallfit", "variable")datat$value<-ifelse(datat$variable=="Depth (m)",datat$depth,datat$time11)# turn the variable depth into the real depth
datat$value<-ifelse(datat$variable=="Time",datat$value/(60*60*24),datat$value)# turn the variable time into the real time in days
datat$variable1<-ifelse(datat$variable=="Time","Time since regular echolocation clicks were resumed (days)",datat$variable)# new labels
datat$whale<-factor(datat$whale,levels=c("Freya","Thora","Mára","Frida","Eistla","Balder","Jonas","Mutti"))# correct order of the whale id2.1.4 Plotting model results/ Figure 5
ggplot(data=datat,aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+facet_wrap(~variable1,scales="free_x", labeller = labeller(variable1 = label_wrap_gen(32)))+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
p1<-ggplot(data=subset(datat,variable1=="Time since regular echolocation clicks were resumed (days)"),aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("")+xlab("Time since regular echolocation clicks \nwere resumed (days)")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
theme(axis.text.y=element_blank(),
# axis.ticks.y=element_blank(),
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
# axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))+ylim(2,12)
p2<-ggplot(data=subset(datat,variable1=="Depth (m)"),aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))+ylim(2,12)
ggpubr::ggarrange(p2,p1,ncol=2,common.legend = T,widths = c(0.5,0.4), legend = "right")
2.1.5 Predicting number of clicks for the clicking periods from the original dataset
onedata$prediction1<-predict(gamwhales1,newdata=(onedata),type="response")# get the prediction from the model2.1.6 Plot with the predicted number of clicks for the clicking periods from the original dataset / Figure 6
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction1),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 17,angle=60, hjust=1),
axis.text.y = element_text(size = 17),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=18),
legend.text= element_text(size=17),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction1),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 17,angle=60, hjust=1),
axis.text.y = element_text(size = 17),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=18),
legend.text= element_text(size=17),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
2.2 GAM number of clicks per second as a function of depth by whale and as a funtion of time by whale
gamwhales11<-gam(clicks~s(depth,k=7)+s(time11,k=1),data=sampledata,family=poisson(link="log"))## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible2.2.1 Summary of the model
summary(gamwhales1)##
## Family: poisson
## Link function: log
##
## Formula:
## clicks ~ s(depth, k = 7, by = as.factor(whale)) + s(time11, k = 1,
## by = as.factor(whale))
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.39154 0.02194 63.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(depth):as.factor(whale)Freya 4.325 4.940 17.035 0.003860 **
## s(depth):as.factor(whale)Thora 5.644 5.940 21.738 0.000932 ***
## s(depth):as.factor(whale)Mára 2.210 2.671 1.631 0.598642
## s(depth):as.factor(whale)Frida 1.000 1.000 19.735 9.49e-06 ***
## s(depth):as.factor(whale)Eistla 2.349 2.919 16.056 0.000969 ***
## s(depth):as.factor(whale)Balder 4.242 4.999 33.387 3.60e-06 ***
## s(depth):as.factor(whale)Jonas 4.750 5.447 29.933 4.59e-05 ***
## s(depth):as.factor(whale)Mutti 1.000 1.000 21.955 3.52e-06 ***
## s(time11):as.factor(whale)Freya 1.729 1.925 13.614 0.005081 **
## s(time11):as.factor(whale)Thora 1.730 1.926 4.325 0.075160 .
## s(time11):as.factor(whale)Mára 1.000 1.000 0.399 0.527811
## s(time11):as.factor(whale)Frida 1.396 1.634 24.510 9.76e-05 ***
## s(time11):as.factor(whale)Eistla 1.000 1.000 2.812 0.093580 .
## s(time11):as.factor(whale)Balder 1.000 1.000 17.759 2.56e-05 ***
## s(time11):as.factor(whale)Jonas 1.001 1.002 0.159 0.690277
## s(time11):as.factor(whale)Mutti 1.000 1.000 0.369 0.543610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0464 Deviance explained = 6.19%
## UBRE = 1.8611 Scale est. = 1 n = 12032.2.2 AIC
AIC(gamwhales1)## [1] 6732.8582.2.3 Predicting influence of the different expl. variables
dd1<-sampledata%>%
summarise(depth=mean(depth,na.rm=T))# creating a data set to predict the influence of time
dd2<-sampledata%>%
summarise(time11=mean(time11,na.rm=T))# creating a data set to predict the influence of depth
originaldata011<-data.frame(depth=sampledata$depth,time11=mean(sampledata$time11,na.rm=T))# creating a data set to predict the influence of depth
originaldata111<-data.frame(time11=sampledata$time11,depth=mean(sampledata$depth,na.rm=T))# creating a data set to predict the influence of time
datat011<-data.frame(pclick=predict(gamwhales11,newdata = originaldata011,type="response"),depth=originaldata011$depth)# predicting the influence of depth in the number of clicks produced
confintervals<-data.frame(confintervals=predict(gamwhales11,newdata=(originaldata011),type="response",se.fit=T))# get the corresponding CI
datat011$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat011$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit
datat011$variable<-"Depth (m)"# which variable corresponds
datat111<-data.frame(pclick=predict(gamwhales11,newdata = originaldata111,type="response"),time11=originaldata111$time11)# predicting the influence of time in the number of clicks produced
confintervals1<-data.frame(confintervals=predict(gamwhales11,newdata=(originaldata111),type="response",se.fit=T))# get the corresponding CI
datat111$bigfit<-confintervals1$confintervals.fit+confintervals1$confintervals.se.fit
datat111$smallfit<-confintervals1$confintervals.fit-confintervals1$confintervals.se.fit
datat111$variable<-"Time"# which variable corresponds
datat01<-full_join(datat011,datat111)#Joinning both infos about the variables ## Joining, by = c("pclick", "bigfit", "smallfit", "variable")datat01$value<-ifelse(datat01$variable=="Depth (m)",datat01$depth,datat01$time11)# turn the variable depth into the real depth
datat01$value<-ifelse(datat01$variable=="Time",datat01$value/(60*60*24),datat01$value)# turn the variable time into the real time in days
datat01$variable1<-ifelse(datat01$variable=="Time","Time since regular echolocation clicks were resumed (days)",datat01$variable)# new labels2.2.4 Plotting model results/ Figure 5
ggplot(data=datat01,aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+facet_wrap(~variable1,scales="free_x", labeller = labeller(variable1 = label_wrap_gen(32)))+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
p11<-ggplot(data=subset(datat01,variable1=="Time since regular echolocation clicks were resumed (days)"),aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("")+xlab("Time since regular echolocation clicks \nwere resumed (days)")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))+ylim(2,8)
p21<-ggplot(data=subset(datat01,variable=="Depth (m)"),aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))+ylim(2,8)
ggpubr::ggarrange(p21,p11,ncol=2,common.legend = T,widths = c(0.5,0.45))
2.2.5 Predicting number of clicks for the clicking periods from the original dataset
onedata$prediction11<-predict(gamwhales11,newdata=(onedata),type="response")# get the prediction from the model2.2.6 Plot with the predicted number of clicks for the clicking periods from the original dataset
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction11),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 17,angle=60, hjust=1),
axis.text.y = element_text(size = 17),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=18),
legend.text= element_text(size=18),
plot.title = element_text(size=21),
strip.text.x = element_text(size = 21))
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction11),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 17,angle=60, hjust=1),
axis.text.y = element_text(size = 17),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=18),
legend.text= element_text(size=18),
plot.title = element_text(size=21),
strip.text.x = element_text(size = 21))
2.3 GAM number of clicks per second as a function of depth by whale
gamwhales12<-gam(clicks~s(depth,k=7,by=as.factor(whale)),data=sampledata,family=poisson(link="log"))2.3.1 Summary of the model
summary(gamwhales12)##
## Family: poisson
## Link function: log
##
## Formula:
## clicks ~ s(depth, k = 7, by = as.factor(whale))
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.41232 0.01765 80.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(depth):as.factor(whale)Freya 5.421 5.814 18.145 0.003678 **
## s(depth):as.factor(whale)Thora 5.764 5.973 25.529 0.000182 ***
## s(depth):as.factor(whale)Mára 2.317 2.867 4.066 0.331349
## s(depth):as.factor(whale)Frida 2.428 2.959 19.767 0.000158 ***
## s(depth):as.factor(whale)Eistla 2.149 2.683 12.328 0.004252 **
## s(depth):as.factor(whale)Balder 4.327 5.078 31.131 1.07e-05 ***
## s(depth):as.factor(whale)Jonas 4.850 5.523 28.508 9.85e-05 ***
## s(depth):as.factor(whale)Mutti 1.000 1.000 19.633 9.96e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0335 Deviance explained = 4.71%
## UBRE = 1.8933 Scale est. = 1 n = 12032.3.2 AIC
AIC(gamwhales12)## [1] 6771.6212.3.3 Predicting influence of the different expl. variables
originaldata112<-data.frame(depth=sampledata$depth,whale=sampledata$whale)# creating a data set to predict the influence of depth
datat12<-data.frame(pclick=predict(gamwhales12,newdata = originaldata112,type="response"),depth=originaldata112$depth,whale=originaldata112$whale)# predicting the influence of depth in the number of clicks produced
confintervals<-data.frame(confintervals=predict(gamwhales12,newdata=(originaldata112),type="response",se.fit=T))# get the corresponding CI
datat12$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat12$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit
datat12$variable<-"Depth (m)"# which variable corresponds
datat2<-(datat12)#Joinning both infos about the variables
datat2$value<-ifelse(datat2$variable=="Depth (m)",datat2$depth,datat2$time11)# turn the variable depth into the real depth
datat2$whale<-factor(datat2$whale,levels=c("Freya","Thora","Mára","Frida","Eistla","Balder","Jonas","Mutti"))# correct order of the whale id2.3.4 Plotting model results
ggplot(data=subset(datat2,variable=="Depth (m)"),aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
2.3.5 Predicting number of clicks for the clicking periods from the original dataset
onedata$prediction12<-predict(gamwhales12,newdata=(onedata),type="response")# get the prediction from the model2.3.6 Plot with the predicted number of clicks for the clicking periods from the original dataset
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction12),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 17),
axis.text.x = element_text(size = 18,angle=60, hjust=1),
axis.text.y = element_text(size = 18),
axis.title.y = element_text(size = 17),
legend.title = element_text(size=17),
legend.text= element_text(size=17),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction12),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 18,angle=60, hjust=1),
axis.text.y = element_text(size = 18),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=17),
legend.text= element_text(size=17),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
2.4 GAM number of clicks per second as a function of depth
gamwhales13<-gam(clicks~s(depth,k=7),data=sampledata,family=poisson(link="log"))2.4.1 Summary of the model
summary(gamwhales13)##
## Family: poisson
## Link function: log
##
## Formula:
## clicks ~ s(depth, k = 7)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.44647 0.01405 103 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(depth) 3.667 4.466 78.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0237 Deviance explained = 2.33%
## UBRE = 1.9234 Scale est. = 1 n = 12032.4.2 AIC
AIC(gamwhales13)## [1] 6807.832.4.3 Predicting influence of the different expl. variables
originaldata113<-data.frame(depth=sampledata$depth)# creating a data set to predict the influence of depth
datat13<-data.frame(pclick=predict(gamwhales13,newdata = originaldata113,type="response"),depth=originaldata113$depth)# predicting the influence of depth in the number of clicks produced
confintervals<-data.frame(confintervals=predict(gamwhales13,newdata=(originaldata113),type="response",se.fit=T))# get the corresponding CI
datat13$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat13$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit
datat13$variable<-"Depth (m)"# which variable corresponds
datat3<-(datat13)#Joinning both infos about the variables
datat3$value<-ifelse(datat3$variable=="Depth (m)",datat3$depth,datat3$time11)# turn the variable depth into the real depth2.4.4 Plotting model results
ggplot(data=subset(datat3,variable=="Depth (m)"),aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
theme_bw()+ scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+ scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
theme(
axis.title.x = element_text(size = 20),
axis.text.x = element_text(size = 19),
axis.text.y = element_text(size = 19),
axis.title.y = element_text(size = 20),
legend.title = element_text(size=20),
legend.text= element_text(size=19),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
2.4.5 Predicting number of clicks for the clicking periods from the original dataset
onedata$prediction13<-predict(gamwhales13,newdata=(onedata),type="response")# get the prediction from the model2.4.6 Plot with the predicted number of clicks for the clicking periods from the original dataset
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction13),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 17,angle=60, hjust=1),
axis.text.y = element_text(size = 17),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=18),
legend.text= element_text(size=17),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
geom_point(color="black",shape=".") +
geom_point(data=onedata,aes(colour=(prediction13),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
na.value = "black")+
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 17,angle=60, hjust=1),
axis.text.y = element_text(size = 17),
axis.title.y = element_text(size = 18),
legend.title = element_text(size=18),
legend.text= element_text(size=17),
plot.title = element_text(size=20),
strip.text.x = element_text(size = 21))
2.5 Comparing models
data.frame("Response variable"=rep("Number of clicks per second",4),"Explanatory variables"=c("depth by narwhal ID and time by narwhal","depth and time","depth by narwhal ID","depth"),"AIC"=c(round(AIC(gamwhales1),3),round(AIC(gamwhales11),3),round(AIC(gamwhales12),3),round(AIC(gamwhales13),3)))2.5.1 Estimating the cue rate using the model results
tagdur<-originaldata%>%
dplyr::group_by(whale)%>%
dplyr::summarise(sumone=unique(len_click),ll=unique(tot_len))
onedata<-full_join(onedata,tagdur)## Joining, by = "whale"est_cr<-onedata%>%
group_by(whale)%>%
summarise(meany=mean(prediction1,na.rm=T),sumone=unique(sumone),ll=unique(ll),mean=(meany*sumone)/ll,sd=sd(prediction1,na.rm=T),cv=sd/meany,
lower = meany-qt(0.975,df=n()-1)*sd/sqrt(n()),
upper = meany+qt(0.975,df=n()-1)*sd/sqrt(n()))2.5.1.1 Results
n <- as.numeric(nrow(est_cr))
xbar <- mean(est_cr$mean,na.rm=T)
s <- sd(est_cr$mean,na.rm=T)
#calculate margin of error
margin <- qt(0.975,df=n-1)*s/sqrt(n)
#calculate lower and upper bounds of confidence interval
low <- xbar - margin
high <- xbar + margin
cv0<-s/xbar
data0<-data.frame(mean=xbar,low=low,high=high,cv=cv0,sd=s,id="Standard Average")
#var(df5$mean)The standard average for the cue rate is 1.412 and the corresponding 95% CI is [0.998 ,1.827]. The cv obtained is 0.351 and the sd is 0.495.
2.5.1.2 Weighted Results
wxbar<- weighted.mean(est_cr$mean,est_cr$ll,na.rm=T)
ws <- sqrt(var.wtd.mean.cochran(est_cr$mean,est_cr$ll))
#calculate margin of error
wmargin <- qt(0.975,df=n-1)*ws
#calculate lower and upper bounds of confidence interval
wlow <- wxbar - wmargin
whigh <- wxbar + wmargin
wcv<-ws/wxbar
#var(df5$mean)
wdata0<-data.frame(mean=wxbar,low=wlow,high=whigh,cv=wcv,sd=ws,id="Weighted Average")The weighted average for the cue rate is 1.29 and the corresponding 95% CI is [0.961 ,1.62]. The cv obtained is 0.108 and the weighted sd is 0.139.
2.5.1.3 Joining results from cue rates
adata<-full_join(data0,wdata0)## Joining, by = c("mean", "low", "high", "cv", "sd", "id")adata$pred<-"Prediction"
data$pred<-"Original"
cr_data<-full_join(adata,data)## Joining, by = c("mean", "low", "high", "cv", "sd", "id", "pred")cr_data<-data.frame(cr_data[cr_data$id%in%c("Weighted Average","Standard Average"),])data.frame(cr_data[cr_data$id%in%c("Weighted Average","Standard Average"),c("id","pred","mean","sd","cv","low","high")])2.5.1.4 Plotting results
ggplot(cr_data, aes(x= id, y= mean,color=pred,fill=pred))+
# geom_point()+
geom_pointrange(aes(ymin=low,ymax=high,color=pred,fill=pred),
position = position_dodge(width = 0.5))+
ylab("Mean Cue Rate")+
xlab("")+
labs(color="",fill="")+
theme(axis.text.x=element_text(angle=30, hjust=1))
2.5.2 Saving the best model and its results
save(file="model_results_n_clicks.rda",list=c("gamwhales1","datat"))